USE ElliDef
USE matrix_oper
USE NodeInfoDef
USE metric
USE Elliptic_Solvers
IMPLICIT NONE
TYPE (NodeInfo) :: Info


INTEGER error,N,nmax
REAL (KIND=qPrec) :: xUpper(1:3),dx,dy,dz,xmax
REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:) :: p,derp,u,v,w,div,pl,dummy,dummy1
REAL (KIND=qPrec), ALLOCATABLE, DIMENSION (:,:) :: ut,vt,wt


INTEGER i,j,k,nx,ny,nz,l
INTEGER, PARAMETER :: solver=PARDISO_SOLVER
!INTEGER, PARAMETER :: solver=IT_FOM
LOGICAL iterate









! populate Info structures

ALLOCATE(Info%ElliInfo)
xmax=25.;nmax=50
Info%mX=(/nmax,10,1/)
nx=Info%mX(1);ny=Info%mX(2);nz=Info%mX(3)
Info%nDim=2
Info%ElliInfo%ndim=Info%ndim
Info%ElliInfo%mX=Info%mX ! last dim has to be one
!!$ALLOCATE(Info%ElliInfo%p(0:Info%mX(1)+1,0:Info%mX(2)+1,0:Info%mX(3)+1))
!!$ALLOCATE(Info%ElliInfo%pold(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3)))
!!$ALLOCATE(Info%ElliInfo%qe(nx+1,ny+1,nz+1,2))
!!$ALLOCATE(Info%ElliInfo%div(nx,ny,nz))

!xUpper(1:3)=(/20.d0,30.d0,1.d0/)

!!$Info%xLower=(/-5000./2200.,0.,0./)
!!$
!!$Info%xUpper=(/15000./2200,1.,1./)

Info%xLower=(/-xmax,0.d0,0.d0/)
!!$
Info%xUpper=(/xmax,1.d0,1.d0/)


Info%dX(1:2)=(Info%xUpper(1:2)-Info%xLower(1:2))/Info%mX(1:2);Info%dX(3)=1
Info%ElliInfo%dX(1:2)=Info%dX(1:2)
Info%ElliInfo%dX(3)=Info%dX(3) ! just to have it set to something
Info%dmax=0.

!Info%mthbc=InternalBC ! assume all Neumann conditions
Info%ElliInfo%mthbc(1:6)=0
!Info%mthbc(5)=InternalBC




CALL setprob(Info)
PRINT*,'Dimension of problem',Info%mX

ALLOCATE(ut(1:ny,1:nz),vt(1:nx,1:nz),wt(0:nx+1,0:ny+1))
CALL DEFINE_TOPO(Info) ! define the quantities used to calculate the metric
error=CREATE_OPERATORS_THIN(Info%ElliInfo)

! zero BC here to be safe 
DO i=1,2*Info%nDim
   Info%ElliInfo%bc(i)%p=zero
END DO

IF(.NOT.check(Info%ElliInfo,CHECK_BC)) THEN
   PRINT*,'BC are not consistent'
   STOP
END IF






! calculates divergence


u=>Info%ElliInfo%qe(1:nx+1,1:ny,1:nz,1);
v=>Info%ElliInfo%qe(1:nx,1:ny+1,1:nz,2);


!!$CALL RANDOM_NUMBER(u)
!!$CAlL RANDOM_NUMBER(v)
u=zero;v=zero
v(:,:,1)=SPREAD(-one*.5*(cos(8.*atan(1.d0)*((/1:nx/)-0.5)/float(nx))),DIM=2,NCOPIES=ny+1)
ut=zero;vt=zero
!!$CALL RANDOM_NUMBER(ut)
!!$CALL RANDOM_NUMBER(vt)



        



! only need OP= optional to use IT wth PARDISO precon
i=POISSON_SOLVERS_THIN(Info%ElliInfo,INIT)!,OP=CSR_PRECOND)
i=POISSON_SOLVERS_THIN(Info%ElliInfo,INIT,OP=8)



ITERATE=.true.;i=1
DO WHILE(ITERATE)
PRINT*,'iteration ',i

!!$
wt=zero
! ********************* 1st vertical **********************
5 Info%ElliInfo%bc(1)%p(1:ny,1:nz)=u(1,:,:)-SUM(u(1,:,1),1)/ny-(ut-SUM(ut(:,1),1)/ny) ! u*-utarget
Info%ElliInfo%bc(2)%p(1:ny,1:nz)=u(nx+1,:,:)-SUM(u(1,:,1),1)/ny-(ut-SUM(ut(:,1),1)/ny)
!!$Info%ElliInfo%bc(3)%p(1:nx,1:nz)=zero
!!$Info%ElliInfo%bc(4)%p(1:nx,1:nz)=zero
Info%ElliInfo%bc(3)%p(1:nx,1:nz)=v(1:nx,1,1:nz)-vt
   Info%ElliInfo%bc(4)%p(1:nx,1:nz)=v(1:nx,ny+1,1:nz)-vt

CALL take_div(Info%ElliInfo,.true.)
div=>Info%ElliInfo%div



! modify upper velocity and upper boundary condition
Info%ElliInfo%bc(4)%p(1:nx,1:nz)=Info%ElliInfo%bc(4)%p(1:nx,1:nz)+ &
     SUM(div,2)/ny
v(:,ny+1,1)=v(:,ny+1,1)-SUM(div(:,:,1),2)/ny
div=div-SPREAD(SUM(div,2)/ny,DIM=2,NCOPIES=ny)
PRINT*,'divergence stat'
!!$
PRINT*,'max divergence before = ',maxval(abs(div))
PRINT*,'rms divergence before = ',sqrt(SUM(div**2)/(nx*ny*nz))
PRINT*,'maxval vert sum div = ',maxval(SUM(div,2))
Info%ElliInfo%p=zero
! solve 1d problems
j=POISSON_SOlVERS_THIN(Info%ElliInfo,SOLVE,OP=8)

! the following is necessary because we need to add off diagonal boundary
! conditions (depending on solution at 0th level) to the 1st order correction



Info%ElliInfo%p(1:nx,1:ny,1)=Info%ElliInfo%p(1:nx,1:ny,1)+wt(1:nx,1:ny)!SPREAD(wt(1:nx,1),DIM=2,NCOPIES=ny)

dummy=>adv(Info%ElliInfo%G(0)%crsptr,Info%ElliInfo%p(1:nx,1:ny,1:nz),(/2,2,2/))

 Info%ElliInfo%p(0:nx+1,0:ny+1,0:nz+1)=Info%ElliInfo%p(0:nx+1,0:ny+1,0:nz+1)+dummy(1:nx+2,1:ny+2,1:nz+2)
DEALLOCATE(dummy)
Info%ElliInfo%p(1:nx,1:ny,1)=Info%ElliInfo%p(1:nx,1:ny,1)-wt(1:nx,1:ny)!SPREAD(wt(1:nx,1),DIM=2,NCOPIES=ny)
! only the off diagonal are needed...
Info%ElliInfo%p(0,:,1)=Info%ElliInfo%p(0,:,1)-wt(1,1:ny)
Info%ElliInfo%p(nx+1,:,1)=Info%ElliInfo%p(nx+1,:,1)-wt(nx,1:ny)
Info%ElliInfo%p(:,0,1)=Info%ElliInfo%p(:,0,1)-wt(:,1)
Info%ElliInfo%p(:,ny+1,1)=Info%ElliInfo%p(:,ny+1,1)-wt(:,ny)
Info%ElliInfo%p(:,:,1)=Info%ElliInfo%p(:,:,1)-SPREAD(Info%ElliInfo%p(:,1,1),DIM=2,NCOPIES=ny+2)
CALL project(Info%ElliInfo,.false.);
wt=Info%ElliInfo%p(:,:,1)

! ************************** end 1th vertical ******************
!
! *************************** 1th horizontal *******************
! BC are now enforced on the average
Info%ElliInfo%p=0
DO j=1,2
   Info%ElliInfo%bc(j)%p(:,1)=SUM(u(1+j/2*nx,:,1),1)/ny-SUM(ut(:,1),1)/ny
ENDDO
   Info%ElliInfo%bc(3)%p(1:nx,1:nz)=v(1:nx,1,1:nz)-vt
   Info%ElliInfo%bc(4)%p(1:nx,1:nz)=v(1:nx,ny+1,1:nz)-vt

PRINT*,'calculates divergence'

CALL take_div(Info%ElliInfo,.true.)
div=>Info%ElliInfo%div

PRINT*,'divergence stat'
!!$
PRINT*,'max divergence before = ',maxval(abs(div))
PRINT*,'rms divergence before = ',sqrt(SUM(div**2)/(nx*ny*nz))
! average vertically divergence
div(:,1,1)=SUM(div(:,:,1),2)/float(ny)
! solve D-1 problem
j=POISSON_SOLVERS_THIN(Info%ElliInfo,SOLVE)
! project the field using this pressure
! first add boundary conditions horizontally
dummy=>adv(Info%ELLIINFO%G(10)%crsptr,Info%ElliInfo%p(1:nx,1:1,1:1),(/2,2,2/))

Info%ElliInfo%p(:,1,1)=Info%ElliInfo%p(:,1,1) +dummy(:,1,1)

DEALLOCATE(dummy);NULLIFY(dummy)

Info%ElliInfo%p(:,:,1)=SPREAD(Info%ElliInfo%p(:,1,1),DIM=2,NCOPIES=ny+2)

CALL project(Info%ElliInfo,.false.);
wt=wt+Info%ElliInfo%p(:,:,1)
! ********************************* end 1th horizontal ***************

CALL take_div(Info%ElliInfo,.false.)
ITERATE=sqrt(SUM(Info%ElliInfo%div**2)/(nx*ny*nz))>1e-8.AND.i<10
IF(ITERATE) THEN

   i=i+1
   PRINT*,'at end of iteration',i-1, 'l2 norm of div is ',sqrt(SUM(Info%ElliInfo%div**2)/(nx*ny*nz))
   GOTO 5
ENDIF

END DO

GOTO 10
div=>adv(Info%ElliInfo%G(4)%crsptr,u,(/-1,0,0/))
v(:,1,1)=vt(:,1)
DO j=2,ny+1
   v(:,j,1)=v(:,j-1,1)-div(:,j-1,1)*Info%dX(2)
ENDDO


PRINT*,v(1:2,:,1)

i=POISSON_SOLVERS_THIN(Info%ElliInfo,FINISH)
Info%ElliInfo%p(:,:,1)=wt

PRINT*,'parison_mopup'


10 open(10,file='u.dat',form='unformatted')
write(10) Info%ElliInfo%qe(1:nx+1,1:ny,1:nz,1)
close(10)
open(10,file='w.dat',form='unformatted')
write(10) Info%ElliInfo%qe(1:nx,1:ny+1,1:nz,2)
close(10)
open(10,file='p.dat',form='unformatted')
write(10) wt(0:nx+1,0:ny+1)
close(10)
open(10,file='div.dat',form='unformatted')
write(10) Info%ElliInfo%div(1:nx,1:ny,1:nz)
close(10)

error=FREEUP_OPERATORS(Info%ElliInfo)
 
DEALLOCATE(Info%ElliInfo)
DEALLOCATE(Info%topo%h,Info%topo%hx)
DEALLOCATE(Info%topo)
DEALLOCATE(ut,vt)


END

